The starting point for hydrograph analysis is to obtain the source data. Let’s see how it goes with sample spas-zagorye.txt discharge data for \(1956-2020\) year range provided with grwat. This data is for Spas-Zagorye gauge on Protva river in Central European plane:
library(sf) # reading and manipulating spatial data
library(tidyverse) # general data wrangling
library(mapview) # interactive mapping of spatial data
library(ecmwfr) # this is to access ERA5 reanalysis data
library(grwat)
mapviewOptions(fgb = FALSE)
# this is path to sample data installed with grwat
path = system.file("extdata", "spas-zagorye.txt", package = "grwat")
# for your own data just provide the full path:
# path = /this/is/the/path/to/discharge/discharge.csv
hdata = read_delim(path, col_names = c('d', 'm', 'y', 'q'), col_types = 'iiid', delim = ' ') # read gauge data
head(hdata) # see the data
#> # A tibble: 6 x 4
#> d m y q
#> <int> <int> <int> <dbl>
#> 1 1 1 1956 5.18
#> 2 2 1 1956 5.18
#> 3 3 1 1956 5.44
#> 4 4 1 1956 5.44
#> 5 5 1 1956 5.44
#> 6 6 1 1956 5.58This example contains the minimum amount of data needed to use the grwat. The first three columns encode the date, and the last column is the main variable — discharge. If the date is encoded as a single column instead of three separate columns, this is also perfect for grwat as input:
hdata_alt = hdata %>%
transmute(D = lubridate::make_date(y, m, d),
Q = q)
head(hdata_alt)
#> # A tibble: 6 x 2
#> D Q
#> <date> <dbl>
#> 1 1956-01-01 5.18
#> 2 1956-01-02 5.18
#> 3 1956-01-03 5.44
#> 4 1956-01-04 5.44
#> 5 1956-01-05 5.44
#> 6 1956-01-06 5.58A sole discharge data is enough to separate the hydrograph into quickflow and baseflow, but is not sufficient to predict the genesis of quickflow cases. Was it due to rain or thaw? To answer such questions you also need precipitation and temperature data. Ideally, these must be measured at the gauge. But often such data is not available. In this case you need to mine this data from external sources.
One of the ways to obtain the temperature and precipitation data is to use reanalyses such as ERA5. Reanalysis data is arranged as regular grids with specific resolution. In particular, the ERA5 data has \(31\) km or \(0.28125\) degrees resolution. To use such data you must tolerate the fact that none of the reanalysis grid nodes will coincide with your gauge. Instead, you have to use the data, which is either
The last two options are the most adequate. Let’s see how it can be done.
First, we need to read the basin spatial data:
# this is path to sample basin geopackage installed with grwat
path = system.file("extdata", "spas-zagorye.gpkg", package = "grwat")
# for your own data just provide the full path:
# path = /this/is/the/path/to/discharge/basin.shp
basin = st_read(path, layer = 'basin') # read basin region
#> Reading layer `basin' from data source
#> `/private/var/folders/5s/rkxr4m8j24569d_p6nj9ld200000gn/T/RtmpFIW0dE/temp_libpathe4ef26265f1d/grwat/extdata/spas-zagorye.gpkg'
#> using driver `GPKG'
#> Simple feature collection with 1 feature and 7 fields
#> Geometry type: MULTIPOLYGON
#> Dimension: XY
#> Bounding box: xmin: 35.41204 ymin: 54.88195 xmax: 36.84138 ymax: 55.57005
#> Geodetic CRS: WGS 84
gauge = st_read(path, layer = 'gauge') # read gauge point
#> Reading layer `gauge' from data source
#> `/private/var/folders/5s/rkxr4m8j24569d_p6nj9ld200000gn/T/RtmpFIW0dE/temp_libpathe4ef26265f1d/grwat/extdata/spas-zagorye.gpkg'
#> using driver `GPKG'
#> Simple feature collection with 1 feature and 0 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 36.62272 ymin: 55.03669 xmax: 36.62272 ymax: 55.03669
#> Geodetic CRS: WGS 84
mapview(basin) + mapview(gauge)Next, we can buffer the data on the specified distance to catch more reanalysis data. For this we use gr_buffer_geo function, which approximates geographic buffer of the specified radius:
basin_buffer = gr_buffer_geo(basin, 25000)
mapview(basin_buffer, col.regions = 'red') +
mapview(basin)Alternatively you can just buffer the gauge point, though it is less meaningful since you will grab reanalysis data that falls out of gauge’s basin:
gauge_buffer = gr_buffer_geo(gauge, 50000)
mapview(gauge_buffer, col.regions = 'red') +
mapview(gauge)grwat is packaged with daily reanalysis which covers the East European territory of Russia. It has a spatial resolution of \(0.75^{\circ} \times 0.75^{\circ}\), and temporal resolution of \(1\) (one) day. Data sources include:
The coverage of the reanalysis is shown below:

The reanalysis consists of two data files, each file is about 850 Mb in size:

Download this data using the FTP link for using with grwat.
rean = gr_read_rean('/Volumes/Data/Spatial/Reanalysis/grwat/pre_1880-2021.nc',
'/Volumes/Data/Spatial/Reanalysis/grwat/temp_1880-2021.nc') # read reanalysis data
hdata_rean = gr_join_rean(hdata, rean, basin_buffer) # join reanalysis data to hydrological series
#> Joining, by = c("Year", "Month", "Day")
head(hdata_rean$df)
#> Day Month Year Q Temp Prec
#> 1 1 1 1956 5.18 -6.46 0.453
#> 2 2 1 1956 5.18 -11.41 0.825
#> 3 3 1 1956 5.44 -10.74 0.260
#> 4 4 1 1956 5.44 -8.05 0.397
#> 5 5 1 1956 5.44 -11.73 0.102
#> 6 6 1 1956 5.58 -20.13 0.032After reanalysis data are joined you can easily plot a map of the derived spatial configuration:
# plot spatial configuration
m = mapview(basin_buffer, col.regions = 'red') +
mapview(basin) +
mapview(hdata_rean$pts, col.regions = 'black') +
mapview(rean$pts, cex = 1)
box = st_bbox(basin_buffer)
center = st_coordinates(st_centroid(basin_buffer))
#> Warning in st_centroid.sf(basin_buffer): st_centroid assumes attributes are
#> constant over geometries of x
m@map %>% leaflet::setView(center[1], center[2], zoom = 7)To be described
gr_fill_gaps() function allows to fill the gaps in the data using simple linear interpolation between nearest observations. You can limit the maximum gap extent by threshold autocorrelation value (autocorr parameter) or expliсit number of observations (nobserv parameter):
tab = gr_fill_gaps(hdata_rean$df,
nobserv = 10)
#> grwat: filled 0 observations using 10 days window
tab = gr_fill_gaps(hdata_rean$df,
autocorr = 0.7)
#> grwat: filled 0 observations using 5 days windowHow the time window is computed in the second case? When the autocorr parameter is used, the ACF (autocorrelation function) is computed first, and then its values are used to obtain the time shift (in days) during which the autocorrelation is higher or equal to the specified value. For Spas-Zagorye data the time lag for \(ACF \geq 0.7\) is 5 days, whoch can be extracted using the purrr::detect_index() function:
afun = acf(hdata_rean$df$Q)
abline(h = 0.7, col = 'red')
abline(v = purrr::detect_index(afun$acf, ~ .x < 0.7), col = 'blue', lwd = 2)
The basic separation procedure is provided by get_baseflow() function. One of the most commonly used approaches is the method by Lyne-Hollick (1979):
resbase = tab %>%
mutate(Date = lubridate::dmy(paste(Day, Month, Year))) %>%
mutate(Qbase = gr_baseflow(Q, method = 'lynehollick'))
# quick look at the table
head(resbase, 10)
#> # A tibble: 10 x 8
#> Day Month Year Q Temp Prec Date Qbase
#> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <date> <dbl>
#> 1 1 1 1956 5.2 -6.5 0.5 1956-01-01 3.70
#> 2 2 1 1956 5.2 -11.4 0.8 1956-01-02 3.79
#> 3 3 1 1956 5.4 -10.7 0.3 1956-01-03 3.88
#> 4 4 1 1956 5.4 -8.1 0.4 1956-01-04 3.96
#> 5 5 1 1956 5.4 -11.7 0.1 1956-01-05 4.04
#> 6 6 1 1956 5.6 -20.1 0 1956-01-06 4.12
#> 7 7 1 1956 5.6 -16.6 0 1956-01-07 4.19
#> 8 8 1 1956 5.4 -15.2 0 1956-01-08 4.26
#> 9 9 1 1956 5.4 -15.9 0.1 1956-01-09 4.32
#> 10 10 1 1956 5.3 -13.3 0.1 1956-01-10 4.39
resbase %>%
filter(Year == 2020) %>%
ggplot() +
geom_area(aes(Date, Q), fill = 'steelblue') +
geom_area(aes(Date, Qbase), fill = 'orangered')
The separation is mainly parameterized by smoothing parameter which is a = 0.925 by default, and number of passes, which are passes = 3 by default. Changing them affects the shape of a baseflow component:
resbase = tab %>%
mutate(Date = lubridate::dmy(paste(Day, Month, Year))) %>%
mutate(Qbase = gr_baseflow(Q, method = 'lynehollick', a = 0.95, passes = 5))
resbase %>%
filter(Year == 2020) %>%
ggplot() +
geom_area(aes(Date, Q), fill = 'steelblue') +
geom_area(aes(Date, Qbase), fill = 'orangered')
Let’s see how different separation methods act in comparison to each other:
methods = c("maxwell",
"boughton",
"jakeman",
"lynehollick",
"chapman")
plots = lapply(methods, function(m) {
resbase = tab %>%
mutate(Date = lubridate::dmy(paste(Day, Month, Year))) %>%
mutate(Qbase = gr_baseflow(Q, method = m))
resbase %>%
filter(Year == 2020) %>%
ggplot() +
geom_area(aes(Date, Q), fill = 'steelblue') +
geom_area(aes(Date, Qbase), fill = 'orangered') +
labs(title = m)
})
patchwork::wrap_plots(plots, ncol = 2)
Advanced separation aims at revealing the genesis of the quickflow. Was it due to a rain or snowmelting? For this, a joint analysis of discharge, temperature and precipitation time series is performed by specialized algorithm, available from gr_separate() function in grwat package.
The genetic separation of hydrograph is controlled by more than 30 parameters. The names and meaning of these paramters can be learned thanks to gr_help_params() function:
gr_help_params()
#> # A tibble: 31 x 11
#> number name example desc units formula comments pics desc_rus role_1
#> <dbl> <chr> <chr> <chr> <chr> <chr> <chr> <lgl> <chr> <chr>
#> 1 1 mome 11 the fi… <NA> <NA> "not used… NA месяц, с… <NA>
#> 2 2 grad 1.5 rate o… %/day <NA> "if chang… NA интенсив… separ…
#> 3 3 grad1 2 rate o… %/day "\"\"" "\"\"" NA интенсив… <NA>
#> 4 4 kdQgr1 150 maximu… % <NA> "can be d… NA максимал… <NA>
#> 5 5 polmo… 2 the ea… <NA> <NA> "can be d… NA cамый ра… <NA>
#> 6 6 polmo… 5 the la… <NA> <NA> "can be d… NA самый по… <NA>
#> 7 7 polko… 15 amount… <NA> <NA> <NA> NA количест… <NA>
#> 8 8 polko… 25 amount… days <NA> <NA> NA количест… <NA>
#> 9 9 polko… 30 amount… days <NA> "can be d… NA количест… <NA>
#> 10 10 polgr… 10 mean r… % <NA> <NA> NA значения… <NA>
#> # … with 21 more rows, and 1 more variable: role_2 <chr>Since the number of parameters is large (31), they are organized as list, which is then fed into gr_separate() function. First, you get the params using the gr_get_params() function. The only parameter of the function is `reg``, which indicates the region for which the parameters must be extracted. After you got the parameters, they can be changed by accessing the list elements:
# Расчленение
p = gr_get_params(reg = 'Midplain')
p$nPav = 5
p$prodspada = 85Next, you can separate the hydrograph:
sep = gr_separate(tab, p)
head(sep)
#> Date Qin Qbase Quick Qseas Qrain Qthaw Qpb Qtype Temp Prec
#> 1 1956-01-01 5.2 5.2 0 0 0 0 0 0 -6.5 0.5
#> 2 1956-01-02 5.2 5.2 0 0 0 0 0 0 -11.4 0.8
#> 3 1956-01-03 5.4 5.4 0 0 0 0 0 0 -10.7 0.3
#> 4 1956-01-04 5.4 5.4 0 0 0 0 0 0 -8.1 0.4
#> 5 1956-01-05 5.4 5.4 0 0 0 0 0 0 -11.7 0.1
#> 6 1956-01-06 5.6 5.6 0 0 0 0 0 0 -20.1 0.0After the hydrograph is separated, it can be summarized in a set of variables:
vars = gr_summarize(sep)
head(vars)
#> # A tibble: 6 x 57
#> Year Year1 Year2 datestart datepolend PolProd Qy Qmax datemax Qygr
#> <int> <int> <dbl> <date> <date> <int> <dbl> <dbl> <date> <dbl>
#> 1 1956 1956 1957 1956-04-08 1956-05-02 24 18.3 467 1956-04-22 7.83
#> 2 1957 1957 1958 1957-03-31 1957-05-01 31 20.3 460 1957-04-08 8.02
#> 3 1958 1958 1959 1958-04-06 1958-05-06 30 27.4 537 1958-04-21 8.28
#> 4 1959 1959 1960 1959-03-31 1959-04-27 27 27.1 406 1959-04-16 7.60
#> 5 1960 1960 1961 1960-03-30 1960-04-26 27 29.5 406 1960-04-15 9.45
#> 6 1961 1961 1962 1961-03-12 1961-04-14 33 18.8 296 1961-04-10 9.81
#> # … with 47 more variables: Qmmsummer <dbl>, monmmsummer <date>, Qmmwin <dbl>,
#> # nommwin <date>, Q30s <dbl>, date30s1 <date>, date30s2 <date>, Q30w <dbl>,
#> # date30w1 <date>, date30w2 <date>, Q10s <dbl>, date10s1 <date>,
#> # date10s2 <date>, Q10w <dbl>, date10w1 <date>, date10w2 <date>, Q5s <dbl>,
#> # date5s1 <date>, date5s2 <date>, Q5w <dbl>, date5w1 <date>, date5w2 <date>,
#> # Wy <dbl>, Wgr <dbl>, Wpol2 <dbl>, Wpol1 <dbl>, Wpol3 <dbl>, Wpavs2 <dbl>,
#> # Wpavs1 <dbl>, Wpavthaw2 <dbl>, Wpavthaw1 <dbl>, WgrS <dbl>, WS <dbl>,
#> # WgrW <dbl>, WW <dbl>, Qmaxpavs <dbl>, Qmaxpavthaw <dbl>,
#> # datemaxpavs <date>, datemaxpavthaw <date>, SumProd <int>,
#> # DaysPavsSum <int>, WinProd <int>, DaysThawWin <int>, CvWin <dbl>,
#> # CvSum <dbl>, CountPavs <int>, CountThaws <int>These functions from grwat package allow you to:
Graphical functions are based on ggplot2 graphics.
You can plot separations for selected years using gr_plot_sep() function:
gr_plot_sep(sep, 1979) # plot single year
#> Warning: Removed 36 rows containing missing values (position_stack).
gr_plot_sep(sep, c(1994, 2001)) # plot two years sequentially
#> Warning: Removed 24 rows containing missing values (position_stack).

And also multiple years on the same layout:
gr_plot_sep(sep, 2016:2019, # plot four years on the same page
layout = matrix(c(1,2,3,4), ncol=2, byrow = T))
#> Warning in max.default(structure(numeric(0), class = "Date"), na.rm = FALSE): no
#> non-missing arguments to max; returning -Inf
To get the detailed description of available variables you can invoke gr_help_vars():
gr_help_vars()
#> # A tibble: 57 x 19
#> ID Position Width Source Name Units Unitsen Readtype Type Test Desc
#> <dbl> <dbl> <dbl> <dbl> <chr> <chr> <chr> <chr> <chr> <dbl> <chr>
#> 1 1 1 7 1 year_… <NA> <NA> integer inte… 0 Номер …
#> 2 2 8 10 1 Year1 Год Year integer inte… 0 Год, к…
#> 3 3 18 10 1 Year2 Год Year integer inte… 0 Год, к…
#> 4 4 28 15 1 dates… Дата Date Date Date 1 Дата н…
#> 5 5 43 15 1 datep… Дата Date Date Date 1 Дата о…
#> 6 57 0 0 0 PolPr… Дней Days integer inte… 1 Продол…
#> 7 6 58 10 1 Qy м^3/с m^3/s double doub… 1 Средни…
#> 8 7 68 10 1 Qmax м^3/с m^3/s double doub… 1 Максим…
#> 9 8 78 15 1 datem… Дата Date Date Date 1 Дата м…
#> 10 9 93 10 1 Qygr м^3/с m^3/s double doub… 1 Средни…
#> # … with 47 more rows, and 8 more variables: Descen <chr>, Group <chr>,
#> # Winter <dbl>, Chart <chr>, Color <chr>, Order <dbl>, Range <chr>,
#> # Problems <chr>Parameters can be statistically tested using test_variables(df, ..., year = NULL, locale='EN') function. Names of the parameters are passed comma-separated in place of .... They are quoted, so you do not need to pass them as character strings, just write their names:
gr_test_vars(vars, Qmax)
#> Warning: `select_()` was deprecated in dplyr 0.7.0.
#> Please use `select()` instead.
#> $ptt
#> $ptt$Qmax
#>
#> Pettitt's test for single change-point detection
#>
#> data: vl[vl_cmp]
#> U* = 481, p-value = 0.01088
#> alternative hypothesis: two.sided
#> sample estimates:
#> probable change point at time K
#> 15
#>
#>
#>
#> $mkt
#> $mkt$Qmax
#>
#> Mann-Kendall trend test
#>
#> data: vl[vl_cmp]
#> z = -3.3955, n = 64, p-value = 0.0006851
#> alternative hypothesis: true S is not equal to 0
#> sample estimates:
#> S varS tau
#> -587.0000000 29785.0000000 -0.2916775
#>
#>
#>
#> $tst
#> $tst$Qmax
#>
#> Sen's slope
#>
#> data: values[fltr]
#> z = -3.3955, n = 64, p-value = 0.0006851
#> alternative hypothesis: true z is not equal to 0
#> 95 percent confidence interval:
#> -5.418605 -1.560976
#> sample estimates:
#> Sen's slope
#> -3.642183
#>
#>
#>
#> $ts_fit
#> $ts_fit$Qmax
#>
#> Call:
#> mblm::mblm(formula = eval(frml), dataframe = df.theil[fltr, ],
#> repeated = FALSE)
#>
#> Coefficients:
#> (Intercept) Year1
#> 7496.073 -3.642
#>
#>
#>
#> $tt
#> $tt$Qmax
#>
#> Welch Two Sample t-test
#>
#> data: d1 and d2
#> t = 3.6606, df = 19.167, p-value = 0.001643
#> alternative hypothesis: true difference in means is not equal to 0
#> 95 percent confidence interval:
#> 75.52882 276.94261
#> sample estimates:
#> mean of x mean of y
#> 419.7857 243.5500
#>
#>
#>
#> $ft
#> $ft$Qmax
#>
#> F test to compare two variances
#>
#> data: d1 and d2
#> F = 1.2603, num df = 13, denom df = 49, p-value = 0.5375
#> alternative hypothesis: true ratio of variances is not equal to 1
#> 95 percent confidence interval:
#> 0.5776842 3.4623990
#> sample estimates:
#> ratio of variances
#> 1.260324
#>
#>
#>
#> $year
#> Qmax
#> 1970
#>
#> $maxval
#> $maxval$Qmax
#> [1] 780
#>
#>
#> $fixed_year
#> [1] FALSE
#>
#> $pvalues
#> N Variable Change.Year Trend
#> 1 1 Maximum annual discharge during seasonal flood wave 1970 -3.64218
#> M1 M2 MeanRatio sd1 sd2 sdRatio Mann.Kendall Pettitt
#> 1 419.7857 243.55 -42 162.9446 145.144 -10.9 0.00069 0.01088
#> Student Fisher
#> 1 0.00164 0.53749This is an example with three variables selected:
tests = gr_test_vars(vars, Qygr, date10w1, Wpol3)
tests$pvalues
#> N Variable
#> 1 1 Annual groundwater discharge ("baseflow") during water-resources year
#> 2 2 First date of 10-day window discharge during winter
#> 3 3 Seasonal flood runoff (with groundwater and rainwater)
#> Change.Year Trend M1 M2 MeanRatio sd1 sd2 sdRatio
#> 1 1979 0.11219 7.85396 13.30172 69.4 1.84951 2.32226 25.6
#> 2 1989 0.97421 22-Apr 10-Jun 49.0 133.00000 150.00000 12.8
#> 3 1988 -0.10666 10.2846 6.46593 -37.1 5.44084 3.72020 -31.6
#> Mann.Kendall Pettitt Student Fisher
#> 1 0.00000 0.00000 0.00000 0.25596
#> 2 0.47594 0.59462 0.17120 0.49666
#> 3 0.00072 0.00628 0.00206 0.03901If you want to test all parameters, just skip variable names:
tests = gr_test_vars(vars)
tests$year # this is a change year detected for each variable
#> Wy Wpol2 Wgr Wpol1 Wpavs2
#> 1978 1974 1977 1988 1977
#> Wpavs1 Wpavthaw2 Wpavthaw1 WgrW WW
#> 1977 2013 1977 1978 1978
#> WgrS WS Qy datestart Qygr
#> 1977 1977 1978 1970 1979
#> datepolend PolProd Qmax datemax Wpol3
#> 1987 1987 1970 1970 1988
#> Qmmwin nommwin Qmmsummer monmmsummer Q30w
#> 1979 1964 1979 1990 1979
#> date30w1 Q30s date30s1 Q10w date10w1
#> 1977 1979 2000 1979 1989
#> Q10s date10s1 Q5w date5w1 Q5s
#> 1981 2000 1979 1998 1981
#> date5s1 Qmaxpavthaw datemaxpavthaw CountThaws DaysThawWin
#> 2000 2010 1993 1969 1984
#> Qmaxpavs datemaxpavs CountPavs DaysPavsSum CvWin
#> 1960 1965 1980 1970 1983
#> WinProd CvSum SumProd
#> 2000 1991 1988Long-term changes are tested against breaking year, which is calculated for each variable using Pettitt test. However, if you want to use a fixed year, you should pass the desired breaking year into change_year parameter:
tests = gr_test_vars(vars, Qmax, Qygr, change_year = 1987)
tests$ft # Fisher F tests to compare two variances
#> $Qmax
#>
#> F test to compare two variances
#>
#> data: d1 and d2
#> F = 1.2603, num df = 13, denom df = 49, p-value = 0.5375
#> alternative hypothesis: true ratio of variances is not equal to 1
#> 95 percent confidence interval:
#> 0.5776842 3.4623990
#> sample estimates:
#> ratio of variances
#> 1.260324
#>
#>
#> $Qygr
#>
#> F test to compare two variances
#>
#> data: d1 and d2
#> F = 0.6343, num df = 22, denom df = 40, p-value = 0.256
#> alternative hypothesis: true ratio of variances is not equal to 1
#> 95 percent confidence interval:
#> 0.3117119 1.4015893
#> sample estimates:
#> ratio of variances
#> 0.6342958Interannual changes are visualized using gr_plot_vars() function. Its syntax is similar to gr_test_vars() and gr_plot_sep():
gr_plot_vars(vars, Qmax) # plot one selected variable
gr_plot_vars(vars, date10w1, Wpol3) # plot two variables sequentially
#> Warning: Removed 20 rows containing non-finite values (stat_smooth).
#> Warning: Removed 20 rows containing missing values (geom_point).
#> Warning: Removed 2 row(s) containing missing values (geom_path).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).

gr_plot_vars(vars, Qmax, Qygr, date10w1, Wpol3, # plot four variables in matrix layout
layout = matrix(c(1,2,3,4), nrow=2, byrow=TRUE))
#> Warning: Removed 20 rows containing non-finite values (stat_smooth).
#> Warning: Removed 20 rows containing missing values (geom_point).
#> Warning: Removed 2 row(s) containing missing values (geom_path).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
You can add the results of statistical tests to the plot by specifying the result of test_variables() function to the tests parameter. In that case the subtitle with test results will be added, Theil-Sen slope and Pettitt test breaking year are drawn as solid (\(p \leq 0.05\)) or dashed (\(p > 0.05\)) lines:
gr_plot_vars(vars, date10w1, Wpol3, DaysThawWin, Qmaxpavs,
tests = gr_test_vars(vars, date10w1, Wpol3, DaysThawWin, Qmaxpavs)) # add test information
#> Warning: Removed 20 rows containing non-finite values (stat_smooth).
#> Warning: Removed 20 rows containing missing values (geom_point).
#> Warning: Removed 2 row(s) containing missing values (geom_path).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).



Note that
testsparameter ofplot_variables()expects the tests for the same variables as those selected for plotting. If you plot variables A, B, C and supply tests for variables X, Y, Z, they will be added without any warnings, and it is your responsibility to keep them in correspondence with each other
Finally, you can plot all variables by not supplying column names to plot_variables() function. In that case tests (if you want to plot them too) should also be calculated for all variables:
gr_plot_vars(vars, tests = gr_test_vars(vars))Long-term changes are the differences between summarized statistics of one variable calculated for two selected periods. Because these statistics reflect the differences in distributions of parameters, grwat visualizes them as box plots using gr_plot_periods() function. The syntax is similar to gr_plot_vars() except that you must provide either tests or year parameter. If both are supplied then tests is prioritized (you can also supply a fixed year when testing variables:
gr_plot_periods(vars, Qy, year = 1978)
gr_plot_periods(vars, Qy, tests = gr_test_vars(vars, Qy))
Multiple plots can be combined on one page using layout parameter:
gr_plot_periods(vars, Qy, Qmax,
tests = gr_test_vars(vars, Qy, Qmax),
layout = matrix(c(1,2)))
To plot long-term changes for all variables just skip variable names in function call:
gr_plot_periods(vars, tests = gr_test_vars(vars))There is also a small helper function that plots a histogram of minimal discharge month for summer and winter periods:
gr_plot_minmonth(vars, year = 1985)
To render HTML report just pass separation and variables to gr_report() function, and provide the output file name:
report = paste(getwd(), 'Spas-Zagorye.html', sep = '/')
gr_report(sep, vars, output = report)
browseURL(report)
| N | Variable | Change.Year | Trend | M1 | M2 | MeanRatio | sd1 | sd2 | sdRatio | Mann.Kendall | Pettitt | Student | Fisher |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | Annual runoff | 1978 | 0.07966 | 22.72886 | 27.16913 | 19.5 | 7.60958 | 7.52459 | -1.1 | 0.13958 | 0.14553 | 0.03142 | 0.92052 |
| 2 | Seasonal flood runoff (without groundwater) | 1974 | -0.10182 | 10.76172 | 6.14213 | -42.9 | 5.18134 | 3.61374 | -30.3 | NA | 0.00421 | 0.00267 | 0.05912 |
| 3 | Annual groundwater runoff | 1977 | 0.13915 | 8.77161 | 15.83465 | 80.5 | 2.41099 | 5.01904 | 108.2 | 1e-05 | 0 | 0 | 0.00084 |
| 4 | Seasonal flood runoff (with groundwater) | 1988 | -0.10188 | 10.05498 | 6.41846 | -36.2 | 5.15313 | 3.59534 | -30.2 | NA | 0.0061 | 0.00206 | 0.05039 |
| 5 | Rain-flood runoff (without groundwater) | 1969 | 0.00062 | 3.01154 | 2.86082 | -5 | 4.61235 | 3.19159 | -30.8 | NA | 0.75438 | 0.91611 | 0.07746 |
| 6 | Rain-flood runoff (with groundwater) | 1977 | 0.02143 | 5.45264 | 7.32468 | 34.3 | 4.61542 | 5.16170 | 11.8 | NA | 0.223 | 0.15677 | 0.61042 |
| 7 | Thaw-flood runoff (without groundwater) | 2013 | 0.00168 | 1.11322 | 2.55362 | 129.4 | 1.98440 | 2.67068 | 34.6 | 0.67234 | 1 | 0.21075 | 0.22686 |
| 8 | Thaw-flood runoff (with groundwater) | 1977 | 0.02143 | 5.45264 | 7.32468 | 34.3 | 4.61542 | 5.16170 | 11.8 | NA | 0.223 | 0.15677 | 0.61042 |
| 9 | Groundwater runoff during winter | 1978 | 0.0414 | 2.68291 | 5.21297 | 94.3 | 1.51114 | 4.64647 | 207.5 | 0.00021 | 5e-05 | 0.00216 | 0 |
| 10 | Winter low flow runoff | 1978 | 0.04513 | 3.66991 | 6.61922 | 80.4 | 2.96717 | 6.19911 | 108.9 | 0.00218 | 0.00072 | 0.01254 | 0.00062 |
| 11 | Groundwater runoff during summer | 1977 | 0.08197 | 5.10734 | 9.3704 | 83.5 | 1.83627 | 2.94587 | 60.4 | 3e-05 | 0 | 0 | 0.02545 |
| 12 | Summer low flow runoff | 1977 | 0.09643 | 7.58077 | 12.10505 | 59.7 | 4.72252 | 5.61700 | 18.9 | 0.00532 | 0.00294 | 0.00149 | 0.40664 |
| 13 | Annual discharge during water-resources year | 1978 | 0.0141 | 19.64125 | 23.0792 | 17.5 | 6.11827 | 6.52138 | 6.6 | 0.72378 | 0.39157 | 0.04257 | 0.77177 |
| 14 | First date of a seasonal flood wave | 1970 | 365.22694 | 01-Apr | 28-Mar | -4 | 7.00000 | 20.00000 | 185.7 | 0.77616 | 0.59462 | 0.19563 | 0.00047 |
| 15 | Annual groundwater discharge (“baseflow”) during water-resources year | 1979 | 0.11099 | 7.85396 | 13.28411 | 69.1 | 1.84951 | 2.32539 | 25.7 | 0 | 0 | 0 | 0.25314 |
| 16 | Last date of a seasonal flood wave | 1976 | 365.06742 | 03-May | 20-Apr | -13 | 19.00000 | 24.00000 | 26.3 | 0.03894 | 0.07722 | 0.0289 | 0.22037 |
| 17 | Duration of a seasonal flood wave | 1987 | -0.2 | 30.29032 | 22.84848 | -24.6 | 14.86426 | 17.08823 | 15 | 0.00026 | 0.00042 | 0.06739 | 0.44521 |
| 18 | Maximum annual discharge during seasonal flood wave | 1970 | -3.64218 | 419.78571 | 243.49 | -42 | 162.94462 | 145.21444 | -10.9 | 0.00069 | 0.01088 | 0.00164 | 0.53899 |
| 19 | Date of a maximum annual discharge during seasonal flood wave | 1967 | 365.22727 | 16-Apr | 15-Apr | -1 | 5.00000 | 42.00000 | 740 | 0.51975 | 0.55801 | 0.95763 | 0 |
| 20 | Seasonal flood runoff (with groundwater and rainwater) | 1988 | -0.1061 | 10.2846 | 6.4776 | -37 | 5.44084 | 3.75705 | -30.9 | 9e-04 | 0.00689 | 0.00218 | 0.04425 |
| 21 | Minimum monthly discharge during winter | 1979 | 0.1239 | 4.1913 | 9.52927 | 127.4 | 1.98973 | 2.38341 | 19.8 | 0 | 0 | 0 | 0.36874 |
| 22 | Month of a minimum monthly discharge during winter | 1964 | 0 | 12-Jun | 09-May | -34 | 54.00000 | 45.00000 | -16.7 | 0.9386 | 0.71943 | 0.13272 | 0.42549 |
| 23 | Minimum monthly discharge during summer | 1979 | 0.06406 | 5.8 | 9.27561 | 59.9 | 5.62737 | 5.90829 | 5 | 0.00115 | 1e-05 | 0.02419 | 0.82616 |
| 24 | Month of a minimum monthly discharge during summer | 1990 | 0 | 06-Jun | 21-Jun | 15 | 39.00000 | 44.00000 | 12.8 | 0.09557 | 0.14111 | 0.16502 | 0.47698 |
| 25 | Minimum 30-day window discharge during winter | 1979 | 0.1031 | 6.19841 | 10.8255 | 74.6 | 2.29125 | 2.18115 | -4.8 | NA | 0 | 0 | 0.76775 |
| 26 | First date of 30-day window discharge during winter | 1977 | 1.5 | 17-Apr | 14-Jun | 58 | 147.00000 | 163.00000 | 10.9 | 0.0616 | 0.19193 | 0.1694 | 0.64123 |
| 27 | Minimum 30-day window discharge during summer | 1979 | 0.06228 | 6.25773 | 9.75442 | 55.9 | 1.70763 | 2.36964 | 38.8 | NA | 4e-05 | 0 | 0.11114 |
| 28 | First date of 30-day window discharge during summer | 2000 | 365.33333 | 13-Jul | 30-Jul | 17 | 29.00000 | 28.00000 | -3.4 | 0.68392 | 0.24187 | 0.03963 | 0.92991 |
| 29 | Minimum 10-day window discharge during winter | 1979 | 0.10816 | 5.27522 | 10.16122 | 92.6 | 1.73743 | 2.66385 | 53.3 | 0 | 0 | 0 | 0.03518 |
| 30 | First date of 10-day window discharge during winter | 1989 | 0.9881 | 22-Apr | 10-Jun | 49 | 133.00000 | 150.00000 | 12.8 | 0.47594 | 0.59462 | 0.17057 | 0.49787 |
| 31 | Minimum 10-day window discharge during summer | 1981 | 0.06235 | 5.93875 | 8.90868 | 50 | 1.58371 | 1.90763 | 20.5 | NA | 3e-05 | 0 | 0.34881 |
| 32 | First date of 10-day window discharge during summer | 2000 | 365.46154 | 21-Jul | 08-Aug | 18 | 26.00000 | 31.00000 | 19.2 | 0.27383 | 0.07921 | 0.03672 | 0.35453 |
| 33 | Minimum 5-day window discharge during winter | 1979 | 0.11155 | 4.90696 | 9.91122 | 102 | 1.74827 | 2.47677 | 41.7 | 0 | 0 | 0 | 0.08371 |
| 34 | First date of 5-day window discharge during winter | 1998 | 0.68819 | 21-May | 30-Mar | -52 | 146.00000 | 118.00000 | -19.2 | 0.25602 | 0.54022 | 0.12579 | 0.30373 |
| 35 | Minimum 5-day window discharge during summer | 1981 | 0.06276 | 5.81333 | 8.73474 | 50.3 | 1.54071 | 1.86502 | 21 | NA | 2e-05 | 0 | 0.33605 |
| 36 | First date of 5-day window discharge during summer | 2000 | 365.5 | 24-Jul | 11-Aug | 18 | 26.00000 | 30.00000 | 15.4 | 0.21508 | 0.06023 | 0.03352 | 0.38446 |
| 37 | Maximum thaw-flood discharge | 1978 | 0.03541 | 21.84435 | 29.14668 | 33.4 | 38.07410 | 38.38868 | 0.8 | 0.54297 | 1 | 0.47135 | 0.99844 |
| 38 | Date of a maximum thaw-flood discharge | 1993 | 266.27778 | 30-Aug | 04-May | -118 | 146.00000 | 133.00000 | -8.9 | 0.0121 | 0.09772 | 0.00137 | 0.63766 |
| 39 | Number of thaw-flood events | 1969 | -0.09091 | 16.07692 | 11.5098 | -28.4 | 7.46616 | 6.16887 | -17.4 | 0.01896 | 0.13894 | 0.05826 | 0.33884 |
| 40 | Number of days with thaw-flood events | 1984 | -0.53638 | 92.28571 | 74.33333 | -19.5 | 29.40755 | 50.76332 | 72.6 | 0.01493 | 0.008 | 0.08142 | 0.00446 |
| 41 | Maximum rain-flood discharge | 1960 | -0.16821 | 125.91617 | 56.90355 | -54.8 | 98.79488 | 50.38677 | -49 | 0.49053 | 1 | 0.25716 | 0.02792 |
| 42 | Date of a maximum rain-flood discharge | 1965 | 365.9099 | 13-May | 14-Jul | 62 | 35.00000 | 80.00000 | 128.6 | 0.25362 | 0.48353 | 0.00073 | 0.0195 |
| 43 | Number of rain-flood events | 1991 | -0.05 | 12.44828 | 10.2381 | -17.8 | 4.98273 | 5.02896 | 0.9 | 0.24364 | 0.3658 | 0.13095 | 0.94633 |
| 44 | Number of days with rain-flood events | 1970 | 0.05719 | 95.78571 | 107.54 | 12.3 | 38.26060 | 39.61046 | 3.5 | 0.81668 | 1 | 0.32462 | 0.94524 |
| 45 | Relative variation of discharge during winter low flow | 1983 | -0.00039 | 0.37243 | 0.37148 | -0.3 | 0.27190 | 0.31670 | 16.5 | 0.84838 | 1 | 0.98974 | 0.42184 |
| 46 | Duration of winter low flow period | 1968 | -0.22967 | 147.75 | 130.30769 | -11.8 | 60.62272 | 68.16432 | 12.4 | 0.26575 | 0.51707 | 0.39206 | 0.70438 |
| 47 | Relative variation of discharge during summer-autumn low flow | 1991 | -0.0018 | 0.679 | 0.57022 | -16 | 0.32442 | 0.24002 | -26 | 0.35695 | 0.76859 | 0.12886 | 0.10584 |
| 48 | Duration of a summer-autumn low flow period | 1988 | 0.36549 | 197 | 212.8125 | 8 | 47.11756 | 45.36052 | -3.7 | 0.05657 | 0.26671 | 0.17637 | 0.83377 |